﻿using System;
using System.Diagnostics.Contracts;

namespace MathService.DSP.Fourier
{

    [Copyright("alglib", url = "alglib.sources.ru")]
    public static class fft
    {
        /// <summary>
        /// Быстрое комплексное одномерное преобразование Фурье
        /// </summary>
        /// <param name="x">Вектор комплексных значений, преобразуемый в спектр</param>
        /// <remarks>
        /// Размер массива N может быть произвольным числом (составным или простым). 
        /// Составной N будут обработаны с использованием вариативного алгоритма Кули-Тьюки с кешированием.
        /// Массивы, размер которых соответствует малому простому числу преобразуются с помощью жестко сода
        /// (по аналогии с кодом FFTW, но без оптимизации низкого уровня), большое простое число элементов 
        /// обрабатывается с помощью алгоритма Блустейна.
        /// 
        /// Быстрые преобразования для гладких N (только простые множители 2, 3, 5), самый быстрый для степеней 2. 
        /// При N имеющих простые множители большие, чем эти, но порядка меньше, чем N, вычисления 
        /// будут примерно в 4 раза медленнее, чем для близких высоко составных N. 
        /// Когда N является простым, скорость будет в 6 раз меньше.
        /// 
        /// Алгоритмическая сложность O(N*logN) для любых N
        /// </remarks>
        [Copyright("29.05.2009 by Bochkanov Sergey", url = "alglib.sources.ru")]
        public static Complex[] FFT(Complex[] x) //, int n)
        {
            var N = x.Length;
            var plan = new FT_Base.ftplan();

            //
            // Special case: N=1, FFT is just identity transform.
            // After this block we assume that N is strictly greater than 1.
            //
            if(N == 1) return new[] { x[0] };

            //
            // convert input array to the more convinient format
            //
            var buf = new double[2 * N];
            for(var i = 0; i < N; i++)
            {
                buf[2 * i + 0] = x[i].Re;
                buf[2 * i + 1] = x[i].Im;
            }

            //
            // Generate plan and execute it.
            //
            // Plan is a combination of a successive factorizations of N and
            // precomputed data. It is much like a FFTW plan, but is not stored
            // between subroutine calls and is much simpler.
            //
            FT_Base.ftbasegeneratecomplexfftplan(N, plan);
            FT_Base.ftbaseexecuteplan(ref buf, 0, plan);

            //
            // result
            //
            var result = new Complex[N];
            for(var i = 0; i < N; i++)
            {
                result[i].Re = buf[2 * i + 0];
                result[i].Im = buf[2 * i + 1];
            }
            return result;
        }


        /// <summary>
        /// Быстрое обратное комплексное одномерное преобразование Фурье
        /// </summary>
        /// <param name="y">Массив значений спектра</param>
        [Copyright("29.05.2009 by Bochkanov Sergey", url = "alglib.sources.ru")]
        public static Complex[] FFT_Complex_Inverse(Complex[] y) //, int N)
        {
            var N = y.Length;
            //Contract.Assert(apserv.isfinitecvector(a, n), "FFTC1DInv: A contains infinite or NAN values!");

            var result = new Complex[N];
            for(var i = 0; i < N; i++)
                result[i] = new Complex(y[i].Re, -y[i].Im);

            result = FFT(result);//, n);

            for(var i = 0; i < N; i++)
                result[i] = result[i].ComplexConjugate; // / N;
            return result;
        }


        /*************************************************************************
        1-dimensional real FFT.

        Algorithm has O(N*logN) complexity for any N (composite or prime).

        INPUT PARAMETERS
            A   -   array[0..N-1] - real function to be transformed
            N   -   problem size

        OUTPUT PARAMETERS
            F   -   DFT of a input array, array[0..N-1]
                    F[j] = SUM(A[k]*exp(-2*pi*sqrt(-1)*j*k/N), k = 0..N-1)

        NOTE
            F[] satisfies symmetry property F[k] = conj(F[N-k]),  so just one half
        of  array  is  usually needed. But for convinience subroutine returns full
        Complex array (with frequencies above N/2), so its result may be  used  by
        other FFT-related subroutines.


          -- ALGLIB --
             Copyright 01.06.2009 by Bochkanov Sergey
        *************************************************************************/
        /// <summary>
        /// Быстрое одномерное вещественное преобразование Фурье
        /// </summary>
        /// <param name="x">Массив входных значений</param>
        /// <value>Массив комплексных значений спектра</value>
        public static Complex[] FFT(double[] x)
        {
            var N = x.Length;

            //Contract.Assert(N > 0, "FFTR1D: incorrect N!");
            //Contract.Assert(a.Length >= N, "FFTR1D: Length(A)<N!");
            //Contract.Assert(apserv.isfinitevector(a, n), "FFTR1D: A contains infinite or NAN values!");

            var plan = new FT_Base.ftplan();

            //
            // Special cases:
            // * N=1, FFT is just identity transform.
            // * N=2, FFT is simple too
            //
            // After this block we assume that N is strictly greater than 2
            //
            if(N == 1)
                return new[] { new Complex(x[0]) };

            if(N == 2)
                return new[] { new Complex(x[0] + x[1]), new Complex(x[0] - x[1]) };

            //
            // Choose between odd-size and even-size FFTs
            //
            var result = new Complex[N];
            if(N % 2 == 0)
            {
                //
                // even-size real FFT, use reduction to the Complex task
                //
                var n2 = N / 2;
                var buf = new double[N];
                for(var i = 0; i < N; i++) buf[i] = x[i];

                FT_Base.ftbasegeneratecomplexfftplan(n2, plan);
                FT_Base.ftbaseexecuteplan(ref buf, 0, plan);

                var pin = -Consts.pi2 / N;
                var N05 = .5 / N;
                for(var i = 0; i <= n2; i++)
                {
                    var idx = 2 * (i % n2);
                    var hn = new Complex(buf[idx + 0], buf[idx + 1]);
                    idx = 2 * ((n2 - i) % n2);
                    var hmnc = new Complex(buf[idx + 0], -buf[idx + 1]);
                    var v = new Complex(-Math.Sin(pin * i), Math.Cos(pin * i));

                    v = hn + hmnc - v * (hn - hmnc);
                    result[i] = new Complex(v.Re * N05, -v.Im * N05);
                }
                for(var i = n2 + 1; i < N; i++)
                    result[i] = result[N - i].ComplexConjugate;
                return result;
            }
            //
            // use Complex FFT
            //
            for(var i = 0; i < N; i++)
                result[i] = new Complex(x[i]);

            return FFT(result); //, n);
        }


        /*************************************************************************
        1-dimensional real inverse FFT.

        Algorithm has O(N*logN) complexity for any N (composite or prime).

        INPUT PARAMETERS
            F   -   array[0..floor(N/2)] - frequencies from forward real FFT
            N   -   problem size

        OUTPUT PARAMETERS
            A   -   inverse DFT of a input array, array[0..N-1]

        NOTE
            F[] should satisfy symmetry property F[k] = conj(F[N-k]), so just  one
        half of frequencies array is needed - elements from 0 to floor(N/2).  F[0]
        is ALWAYS real. If N is even F[floor(N/2)] is real too. If N is odd,  then
        F[floor(N/2)] has no special properties.

        Relying on properties noted above, FFTR1DInv subroutine uses only elements
        from 0th to floor(N/2)-th. It ignores imaginary part of F[0],  and in case
        N is even it ignores imaginary part of F[floor(N/2)] too.

        When you call this function using full arguments list - "FFTR1DInv(F,N,A)"
        - you can pass either either frequencies array with N elements or  reduced
        array with roughly N/2 elements - subroutine will  successfully  transform
        both.

        If you call this function using reduced arguments list -  "FFTR1DInv(F,A)"
        - you must pass FULL array with N elements (although higher  N/2 are still
        not used) because array size is used to automatically determine FFT length


          -- ALGLIB --
             Copyright 01.06.2009 by Bochkanov Sergey
        *************************************************************************/
        public static void FFT_Real_Inverse(Complex[] f, int n, ref double[] a)
        {

            Contract.Assert(n > 0, "FFTR1DInv: incorrect N!");
            Contract.Assert(f.Length > (int)Math.Floor(n / (double)2), "FFTR1DInv: Length(F)<Floor(N/2)+1!");
            //Contract.Assert(math.isfinite(f[0].x), "FFTR1DInv: F contains infinite or NAN values!");

            //for(var i = 1; i < (int)Math.Floor(n / (double)2); i++)
            //{
            //    Contract.Assert(math.isfinite(f[i].x) && math.isfinite(f[i].y), "FFTR1DInv: F contains infinite or NAN values!");
            //}
            //Contract.Assert(math.isfinite(f[(int)Math.Floor(n / (double)2)].x),
            //          "FFTR1DInv: F contains infinite or NAN values!");
            //if(n % 2 != 0)
            //{
            //    Contract.Assert(math.isfinite(f[(int)Math.Floor(n / (double)2)].y),
            //              "FFTR1DInv: F contains infinite or NAN values!");
            //}

            //
            // Special case: N=1, FFT is just identity transform.
            // After this block we assume that N is strictly greater than 1.
            //
            if(n == 1)
            {
                a = new[] { f[0].Re };
                return;
            }

            //
            // inverse real FFT is reduced to the inverse real FHT,
            // which is reduced to the forward real FHT,
            // which is reduced to the forward real FFT.
            //
            // Don't worry, it is really compact and efficient reduction :)
            //
            var h = new double[n];
            a = new double[n];
            h[0] = f[0].Re;
            for(var i = 1; i < (int)Math.Floor(n / (double)2); i++)
            {
                h[i] = f[i].Re - f[i].Im;
                h[n - i] = f[i].Re + f[i].Im;
            }

            if(n % 2 == 0)
                h[(int)Math.Floor(n / (double)2)] = f[(int)Math.Floor(n / (double)2)].Re;
            else
            {
                h[(int)Math.Floor(n / (double)2)] = f[(int)Math.Floor(n / (double)2)].Re -
                                                    f[(int)Math.Floor(n / (double)2)].Im;
                h[(int)Math.Floor(n / (double)2) + 1] = f[(int)Math.Floor(n / (double)2)].Re +
                                                        f[(int)Math.Floor(n / (double)2)].Im;
            }

            var fh = FFT(h); //, out fh);

            for(var i = 0; i < n; i++)
                a[i] = (fh[i].Re - fh[i].Im) / n;
        }


        /*************************************************************************
        Internal subroutine. Never call it directly!


          -- ALGLIB --
             Copyright 01.06.2009 by Bochkanov Sergey
        *************************************************************************/

        /*
                private static void FFT_Real_InternalEven(ref double[] a, int n, ref double[] buf, FT_Base.ftplan plan)
                {
                    Contract.Assert(n > 0 && n % 2 == 0, "FFTR1DEvenInplace: incorrect N!");

                    //
                    // Special cases:
                    // * N=2
                    //
                    // After this block we assume that N is strictly greater than 2
                    //
                    if(n == 2)
                    {
                        var x = a[0] + a[1];
                        var y = a[0] - a[1];
                        a[0] = x;
                        a[1] = y;
                        return;
                    }

                    //
                    // even-size real FFT, use reduction to the Complex task
                    //
                    var n2 = n / 2;
                    for(var i = 0; i < n; i++)
                        buf[i] = a[i];

                    FT_Base.ftbaseexecuteplan(ref buf, 0, plan);

                    a[0] = buf[0] + buf[1];
                    var pin = -Consts.pi2 / n;
                    for(var i = 1; i < n2; i++)
                    {
                        var idx = 2 * (i % n2);
                        var hn = new Complex(buf[idx + 0], buf[idx + 1]);
                        idx = 2 * (n2 - i);
                        var hmnc = new Complex(buf[idx + 0], -buf[idx + 1]);
                        var v = new Complex(-Math.Sin(pin * i), Math.Cos(pin * i));
                        v = hn + hmnc - v * (hn - hmnc);
                        a[2 * i + 0] = .5 * v.Re;
                        a[2 * i + 1] = .5 * v.Im;
                    }
                    a[1] = buf[0] - buf[1];
                }
        */


        /*************************************************************************
        Internal subroutine. Never call it directly!


          -- ALGLIB --
             Copyright 01.06.2009 by Bochkanov Sergey
        *************************************************************************/

        /*
                private static void FFT_Real_Inverse_InternalEven(ref double[] a, int n, ref double[] buf, FT_Base.ftplan plan)
                {
                    Contract.Assert(n > 0 && n % 2 == 0, "FFTR1DInvInternalEven: incorrect N!");

                    //
                    // Special cases:
                    // * N=2
                    //
                    // After this block we assume that N is strictly greater than 2
                    //
                    if(n == 2)
                    {
                        var x = .5 * (a[0] + a[1]);
                        var y = .5 * (a[0] - a[1]);
                        a[0] = x;
                        a[1] = y;
                        return;
                    }

                    //
                    // inverse real FFT is reduced to the inverse real FHT,
                    // which is reduced to the forward real FHT,
                    // which is reduced to the forward real FFT.
                    //
                    // Don't worry, it is really compact and efficient reduction :)
                    //
                    var n2 = n / 2;
                    buf[0] = a[0];
                    for(var i = 1; i < n2; i++)
                    {
                        var x = a[2 * i + 0];
                        var y = a[2 * i + 1];
                        buf[i] = x - y;
                        buf[n - i] = x + y;
                    }

                    buf[n2] = a[1];
                    FFT_Real_InternalEven(ref buf, n, ref a, plan);

                    a[0] = buf[0] / n;
                    var t = 1.0 / n;
                    for(var i = 1; i < n2; i++)
                    {
                        var x = buf[2 * i + 0];
                        var y = buf[2 * i + 1];
                        a[i] = t * (x - y);
                        a[n - i] = t * (x + y);
                    }
                    a[n2] = buf[1] / n;
                }
        */


        private static class FT_Base
        {
            public const int ftbaseplanentrysize = 8;
            public const int ftbasecffttask = 0;
            public const int ftbaserfhttask = 1;
            public const int ftbaserffttask = 2;
            public const int fftcooleytukeyplan = 0;
            public const int fftbluesteinplan = 1;
            public const int fftcodeletplan = 2;
            public const int fhtcooleytukeyplan = 3;
            public const int fhtcodeletplan = 4;
            public const int fftrealcooleytukeyplan = 5;
            public const int fftemptyplan = 6;
            public const int fhtn2plan = 999;
            public const int ftbaseupdatetw = 4;
            public const int ftbasecodeletrecommended = 5;
            /*
                        public const double ftbaseinefficiencyfactor = 1.3;
            */
            public const int ftbasemaxsmoothfactor = 5;


            /*************************************************************************
            This subroutine generates FFT plan - a decomposition of a N-length FFT to
            the more simpler operations. Plan consists of the root entry and the child
            entries.

            Subroutine parameters:
                N               task size
            
            Output parameters:
                Plan            plan

              -- ALGLIB --
                 Copyright 01.05.2009 by Bochkanov Sergey
            *************************************************************************/
            public static void ftbasegeneratecomplexfftplan(int n, ftplan plan)
            {

                var planarraysize = 1;
                plan.plan = new int[planarraysize];

                var plansize = 0;
                var precomputedsize = 0;
                var stackmemsize = 0;
                var stackptr = 0;
                var tmpmemsize = 2 * n;
                ftbasegenerateplanrec(n, ftbasecffttask, plan, ref plansize, ref precomputedsize, ref planarraysize,
                                      ref tmpmemsize, ref stackmemsize, stackptr);

                Contract.Assert(stackptr == 0, "Internal error in FTBaseGenerateComplexFFTPlan: stack ptr!");

                plan.stackbuf = new double[Math.Max(stackmemsize, 1)];
                plan.tmpbuf = new double[Math.Max(tmpmemsize, 1)];
                plan.precomputed = new double[Math.Max(precomputedsize, 1)];
                stackptr = 0;

                ftbaseprecomputeplanrec(plan, 0, stackptr);
                Contract.Assert(stackptr == 0, "Internal error in FTBaseGenerateComplexFFTPlan: stack ptr!");
            }


            /*************************************************************************
            Generates real FFT plan
            *************************************************************************/
            /*
                        public static void ftbasegeneraterealfftplan(int n, ftplan plan)
                        {
                            var planarraysize = 1;
                            plan.plan = new int[planarraysize];

                            var plansize = 0;
                            var precomputedsize = 0;
                            var stackmemsize = 0;
                            var stackptr = 0;
                            var tmpmemsize = 2 * n;
                            ftbasegenerateplanrec(n, ftbaserffttask, plan, ref plansize, ref precomputedsize, ref planarraysize,
                                                  ref tmpmemsize, ref stackmemsize, stackptr);

                            Contract.Assert(stackptr == 0, "Internal error in FTBaseGenerateRealFFTPlan: stack ptr!");

                            plan.stackbuf = new double[Math.Max(stackmemsize, 1)];
                            plan.tmpbuf = new double[Math.Max(tmpmemsize, 1)];
                            plan.precomputed = new double[Math.Max(precomputedsize, 1)];

                            stackptr = 0;
                            ftbaseprecomputeplanrec(plan, 0, stackptr);
                            Contract.Assert(stackptr == 0, "Internal error in FTBaseGenerateRealFFTPlan: stack ptr!");
                        }
            */


            /*************************************************************************
            Generates real FHT plan
            *************************************************************************/
            /*
                        public static void ftbasegeneraterealfhtplan(int n, ftplan plan)
                        {
                            var planarraysize = 1;
                            plan.plan = new int[planarraysize];

                            var plansize = 0;
                            var precomputedsize = 0;
                            var stackmemsize = 0;
                            var stackptr = 0;
                            var tmpmemsize = n;
                            ftbasegenerateplanrec(n, ftbaserfhttask, plan, ref plansize, ref precomputedsize, ref planarraysize,
                                                  ref tmpmemsize, ref stackmemsize, stackptr);

                            Contract.Assert(stackptr == 0, "Internal error in FTBaseGenerateRealFHTPlan: stack ptr!");

                            plan.stackbuf = new double[Math.Max(stackmemsize, 1)];
                            plan.tmpbuf = new double[Math.Max(tmpmemsize, 1)];
                            plan.precomputed = new double[Math.Max(precomputedsize, 1)];
                            stackptr = 0;

                            ftbaseprecomputeplanrec(plan, 0, stackptr);
                            Contract.Assert(stackptr == 0, "Internal error in FTBaseGenerateRealFHTPlan: stack ptr!");
                        }
            */


            /*************************************************************************
            This subroutine executes FFT/FHT plan.

            If Plan is a:
            * complex FFT plan  -   sizeof(A)=2*N,
                                    A contains interleaved real/imaginary values
            * real FFT plan     -   sizeof(A)=2*N,
                                    A contains real values interleaved with zeros
            * real FHT plan     -   sizeof(A)=2*N,
                                    A contains real values interleaved with zeros

              -- ALGLIB --
                 Copyright 01.05.2009 by Bochkanov Sergey
            *************************************************************************/
            public static void ftbaseexecuteplan(ref double[] a, int aoffset, ftplan plan)
            {
                const int stackptr = 0;
                ftbaseexecuteplanrec(ref a, aoffset, plan, 0, stackptr);
            }


            /*************************************************************************
            Recurrent subroutine for the FTBaseExecutePlan

            Parameters:
                A           FFT'ed array
                AOffset     offset of the FFT'ed part (distance is measured in doubles)

              -- ALGLIB --
                 Copyright 01.05.2009 by Bochkanov Sergey
            *************************************************************************/
            public static void ftbaseexecuteplanrec(ref double[] a, int aoffset, ftplan plan, int entryoffset, int stackptr)
            {
                if(plan.plan[entryoffset + 3] == fftemptyplan)
                    return;

                int i;
                int n1;
                int n2;
                if(plan.plan[entryoffset + 3] == fftcooleytukeyplan)
                {
                    //
                    // Cooley-Tukey plan
                    // * transposition
                    // * row-wise FFT
                    // * twiddle factors:
                    //   - TwBase is a basis twiddle factor for I=1, J=1
                    //   - TwRow is a twiddle factor for a second element in a row (J=1)
                    //   - Tw is a twiddle factor for a current element
                    // * transposition again
                    // * row-wise FFT again
                    //
                    n1 = plan.plan[entryoffset + 1];
                    n2 = plan.plan[entryoffset + 2];

                    internalcomplexlintranspose(ref a, n1, n2, aoffset, ref plan.tmpbuf);

                    for(i = 0; i <= n2 - 1; i++)
                        ftbaseexecuteplanrec(ref a, aoffset + i * n1 * 2, plan, plan.plan[entryoffset + 5], stackptr);

                    ffttwcalc(ref a, aoffset, n1, n2);
                    internalcomplexlintranspose(ref a, n2, n1, aoffset, ref plan.tmpbuf);

                    for(i = 0; i <= n1 - 1; i++)
                        ftbaseexecuteplanrec(ref a, aoffset + i * n2 * 2, plan, plan.plan[entryoffset + 6], stackptr);

                    internalcomplexlintranspose(ref a, n1, n2, aoffset, ref plan.tmpbuf);
                    return;
                }

                int offs;
                int offs2;
                double hk;
                double hnk;
                if(plan.plan[entryoffset + 3] == fftrealcooleytukeyplan)
                {
                    //
                    // Cooley-Tukey plan
                    // * transposition
                    // * row-wise FFT
                    // * twiddle factors:
                    //   - TwBase is a basis twiddle factor for I=1, J=1
                    //   - TwRow is a twiddle factor for a second element in a row (J=1)
                    //   - Tw is a twiddle factor for a current element
                    // * transposition again
                    // * row-wise FFT again
                    //
                    n1 = plan.plan[entryoffset + 1];
                    n2 = plan.plan[entryoffset + 2];

                    internalcomplexlintranspose(ref a, n2, n1, aoffset, ref plan.tmpbuf);


                    for(i = 0; i <= n1 / 2 - 1; i++)
                    {
                        //
                        // pack two adjacent smaller real FFT's together,
                        // make one complex FFT,
                        // unpack result
                        //
                        offs = aoffset + 2 * i * n2 * 2;
                        int k;
                        for(k = 0; k <= n2 - 1; k++)
                            a[offs + 2 * k + 1] = a[offs + 2 * n2 + 2 * k + 0];

                        ftbaseexecuteplanrec(ref a, offs, plan, plan.plan[entryoffset + 6], stackptr);
                        plan.tmpbuf[0] = a[offs + 0];
                        plan.tmpbuf[1] = 0;
                        plan.tmpbuf[2 * n2 + 0] = a[offs + 1];
                        plan.tmpbuf[2 * n2 + 1] = 0;

                        for(k = 1; k <= n2 - 1; k++)
                        {
                            var offs1 = 2 * k;
                            offs2 = 2 * n2 + 2 * k;
                            hk = a[offs + 2 * k + 0];
                            hnk = a[offs + 2 * (n2 - k) + 0];
                            plan.tmpbuf[offs1 + 0] = 0.5 * (hk + hnk);
                            plan.tmpbuf[offs2 + 1] = -(0.5 * (hk - hnk));
                            hk = a[offs + 2 * k + 1];
                            hnk = a[offs + 2 * (n2 - k) + 1];
                            plan.tmpbuf[offs2 + 0] = .5 * (hk + hnk);
                            plan.tmpbuf[offs1 + 1] = .5 * (hk - hnk);
                        }

                        var i1_ = 0 - offs;
                        for(var i_ = offs; i_ <= offs + 2 * n2 * 2 - 1; i_++)
                            a[i_] = plan.tmpbuf[i_ + i1_];
                    }
                    if(n1 % 2 != 0)
                        ftbaseexecuteplanrec(ref a, aoffset + (n1 - 1) * n2 * 2, plan, plan.plan[entryoffset + 6], stackptr);

                    ffttwcalc(ref a, aoffset, n2, n1);
                    internalcomplexlintranspose(ref a, n1, n2, aoffset, ref plan.tmpbuf);

                    for(i = 0; i <= n2 - 1; i++)
                        ftbaseexecuteplanrec(ref a, aoffset + i * n1 * 2, plan, plan.plan[entryoffset + 5], stackptr);

                    internalcomplexlintranspose(ref a, n2, n1, aoffset, ref plan.tmpbuf);
                    return;
                }

                int offsa;
                if(plan.plan[entryoffset + 3] == fhtcooleytukeyplan)
                {
                    //
                    // Cooley-Tukey FHT plan:
                    // * transpose                    \
                    // * smaller FHT's                |
                    // * pre-process                  |
                    // * multiply by twiddle factors  || corresponds to multiplication by H1
                    // * post-process                 |
                    // * transpose again              /
                    // * multiply by H2 (smaller FHT's)
                    // * final transposition
                    //
                    // For more details see Vitezslav Vesely, "Fast algorithms
                    // of Fourier and Hartley transform and their implementation in MATLAB",
                    // page 31.
                    //
                    n1 = plan.plan[entryoffset + 1];
                    n2 = plan.plan[entryoffset + 2];
                    //n = n1 * n2;

                    internalreallintranspose(ref a, n1, n2, aoffset, ref plan.tmpbuf);

                    for(i = 0; i <= n2 - 1; i++)
                        ftbaseexecuteplanrec(ref a, aoffset + i * n1, plan, plan.plan[entryoffset + 5], stackptr);
                    int j;
                    for(i = 0; i <= n2 - 1; i++)
                        for(j = 0; j <= n1 - 1; j++)
                        {
                            offsa = aoffset + i * n1;
                            hk = a[offsa + j];
                            hnk = a[offsa + (n1 - j) % n1];
                            offs = 2 * (i * n1 + j);

                            plan.tmpbuf[offs + 0] = -.5 * (hnk - hk);
                            plan.tmpbuf[offs + 1] = .5 * (hk + hnk);
                        }

                    ffttwcalc(ref plan.tmpbuf, 0, n1, n2);

                    for(j = 0; j <= n1 - 1; j++)
                        a[aoffset + j] = plan.tmpbuf[2 * j + 0] + plan.tmpbuf[2 * j + 1];

                    if(n2 % 2 == 0)
                    {
                        offs = 2 * (n2 / 2) * n1;
                        offsa = aoffset + n2 / 2 * n1;

                        for(j = 0; j <= n1 - 1; j++)
                            a[offsa + j] = plan.tmpbuf[offs + 2 * j + 0] + plan.tmpbuf[offs + 2 * j + 1];
                    }

                    for(i = 1; i <= (n2 + 1) / 2 - 1; i++)
                    {
                        offs = 2 * i * n1;
                        offs2 = 2 * (n2 - i) * n1;
                        offsa = aoffset + i * n1;

                        for(j = 0; j <= n1 - 1; j++)
                            a[offsa + j] = plan.tmpbuf[offs + 2 * j + 1] + plan.tmpbuf[offs2 + 2 * j + 0];

                        offsa = aoffset + (n2 - i) * n1;

                        for(j = 0; j <= n1 - 1; j++)
                            a[offsa + j] = plan.tmpbuf[offs + 2 * j + 0] + plan.tmpbuf[offs2 + 2 * j + 1];
                    }

                    internalreallintranspose(ref a, n2, n1, aoffset, ref plan.tmpbuf);

                    for(i = 0; i <= n1 - 1; i++)
                        ftbaseexecuteplanrec(ref a, aoffset + i * n2, plan, plan.plan[entryoffset + 6], stackptr);

                    internalreallintranspose(ref a, n1, n2, aoffset, ref plan.tmpbuf);
                    return;
                }

                int n;
                if(plan.plan[entryoffset + 3] == fhtn2plan)
                {
                    //
                    // Cooley-Tukey FHT plan
                    //
                    n1 = plan.plan[entryoffset + 1];
                    n2 = plan.plan[entryoffset + 2];
                    n = n1 * n2;

                    reffht(ref a, n, aoffset);
                    return;
                }

                double s1x;
                double s2x;
                double s3y;
                double s4x;
                double s5y;

                double c1;
                double c2;
                double c3;
                double c4;
                double c5;

                double m1x;
                double m2x;
                double m2y;
                double m3y;

                double t1x;
                double t2x;
                double t3x;
                double t4x;
                double t5x;

                double a0x;
                double a1x;
                double a2x;
                double a3x;

                double v0;
                if(plan.plan[entryoffset + 3] == fftcodeletplan)
                {
                    n1 = plan.plan[entryoffset + 1];
                    n2 = plan.plan[entryoffset + 2];
                    n = n1 * n2;

                    double a0y;
                    double a1y;
                    if(n == 2)
                    {
                        a0x = a[aoffset + 0];
                        a0y = a[aoffset + 1];
                        a1x = a[aoffset + 2];
                        a1y = a[aoffset + 3];

                        v0 = a0x + a1x;
                        var v1 = a0y + a1y;
                        var v2 = a0x - a1x;
                        var v3 = a0y - a1y;

                        a[aoffset + 0] = v0;
                        a[aoffset + 1] = v1;
                        a[aoffset + 2] = v2;
                        a[aoffset + 3] = v3;
                        return;
                    }

                    double s1y;
                    double t1y;
                    double a2y;


                    double m1y;
                    if(n == 3)
                    {
                        offs = plan.plan[entryoffset + 7];

                        c1 = plan.precomputed[offs + 0];
                        c2 = plan.precomputed[offs + 1];

                        a0x = a[aoffset + 0];
                        a0y = a[aoffset + 1];
                        a1x = a[aoffset + 2];
                        a1y = a[aoffset + 3];
                        a2x = a[aoffset + 4];
                        a2y = a[aoffset + 5];

                        t1x = a1x + a2x;
                        t1y = a1y + a2y;

                        a0x = a0x + t1x;
                        a0y = a0y + t1y;

                        m1x = c1 * t1x;
                        m1y = c1 * t1y;
                        m2x = c2 * (a1y - a2y);
                        m2y = c2 * (a2x - a1x);

                        s1x = a0x + m1x;
                        s1y = a0y + m1y;

                        a1x = s1x + m2x;
                        a1y = s1y + m2y;
                        a2x = s1x - m2x;
                        a2y = s1y - m2y;

                        a[aoffset + 0] = a0x;
                        a[aoffset + 1] = a0y;
                        a[aoffset + 2] = a1x;
                        a[aoffset + 3] = a1y;
                        a[aoffset + 4] = a2x;
                        a[aoffset + 5] = a2y;
                        return;
                    }

                    double t2y;
                    double m3x;
                    if(n == 4)
                    {
                        a0x = a[aoffset + 0];
                        a0y = a[aoffset + 1];
                        a1x = a[aoffset + 2];
                        a1y = a[aoffset + 3];
                        a2x = a[aoffset + 4];
                        a2y = a[aoffset + 5];
                        a3x = a[aoffset + 6];
                        var a3y = a[aoffset + 7];

                        t1x = a0x + a2x;
                        t1y = a0y + a2y;
                        t2x = a1x + a3x;
                        t2y = a1y + a3y;

                        m2x = a0x - a2x;
                        m2y = a0y - a2y;
                        m3x = a1y - a3y;
                        m3y = a3x - a1x;

                        a[aoffset + 0] = t1x + t2x;
                        a[aoffset + 1] = t1y + t2y;
                        a[aoffset + 4] = t1x - t2x;
                        a[aoffset + 5] = t1y - t2y;
                        a[aoffset + 2] = m2x + m3x;
                        a[aoffset + 3] = m2y + m3y;
                        a[aoffset + 6] = m2x - m3x;
                        a[aoffset + 7] = m2y - m3y;
                        return;
                    }

                    if(n == 5)
                    {
                        offs = plan.plan[entryoffset + 7];

                        c1 = plan.precomputed[offs + 0];
                        c2 = plan.precomputed[offs + 1];
                        c3 = plan.precomputed[offs + 2];
                        c4 = plan.precomputed[offs + 3];
                        c5 = plan.precomputed[offs + 4];

                        t1x = a[aoffset + 2] + a[aoffset + 8];
                        t1y = a[aoffset + 3] + a[aoffset + 9];
                        t2x = a[aoffset + 4] + a[aoffset + 6];
                        t2y = a[aoffset + 5] + a[aoffset + 7];
                        t3x = a[aoffset + 2] - a[aoffset + 8];
                        var t3y = a[aoffset + 3] - a[aoffset + 9];
                        t4x = a[aoffset + 6] - a[aoffset + 4];
                        var t4y = a[aoffset + 7] - a[aoffset + 5];

                        t5x = t1x + t2x;
                        var t5y = t1y + t2y;

                        a[aoffset + 0] = a[aoffset + 0] + t5x;
                        a[aoffset + 1] = a[aoffset + 1] + t5y;

                        m1x = c1 * t5x;
                        m1y = c1 * t5y;
                        m2x = c2 * (t1x - t2x);
                        m2y = c2 * (t1y - t2y);
                        m3x = -(c3 * (t3y + t4y));
                        m3y = c3 * (t3x + t4x);
                        var m4x = -(c4 * t4y);
                        var m4y = c4 * t4x;
                        var m5x = -(c5 * t3y);
                        var m5y = c5 * t3x;

                        var s3x = m3x - m4x;
                        s3y = m3y - m4y;

                        var s5x = m3x + m5x;
                        s5y = m3y + m5y;

                        s1x = a[aoffset + 0] + m1x;
                        s1y = a[aoffset + 1] + m1y;
                        s2x = s1x + m2x;
                        var s2y = s1y + m2y;
                        s4x = s1x - m2x;
                        var s4y = s1y - m2y;

                        a[aoffset + 2] = s2x + s3x;
                        a[aoffset + 3] = s2y + s3y;
                        a[aoffset + 4] = s4x + s5x;
                        a[aoffset + 5] = s4y + s5y;
                        a[aoffset + 6] = s4x - s5x;
                        a[aoffset + 7] = s4y - s5y;
                        a[aoffset + 8] = s2x - s3x;
                        a[aoffset + 9] = s2y - s3y;
                        return;
                    }
                }

                if(plan.plan[entryoffset + 3] == fhtcodeletplan)
                {
                    n1 = plan.plan[entryoffset + 1];
                    n2 = plan.plan[entryoffset + 2];
                    n = n1 * n2;

                    if(n == 2)
                    {
                        a0x = a[aoffset + 0];
                        a1x = a[aoffset + 1];

                        a[aoffset + 0] = a0x + a1x;
                        a[aoffset + 1] = a0x - a1x;
                        return;
                    }

                    if(n == 3)
                    {
                        offs = plan.plan[entryoffset + 7];

                        c1 = plan.precomputed[offs + 0];
                        c2 = plan.precomputed[offs + 1];

                        a0x = a[aoffset + 0];
                        a1x = a[aoffset + 1];
                        a2x = a[aoffset + 2];

                        t1x = a1x + a2x;
                        a0x = a0x + t1x;

                        m1x = c1 * t1x;
                        m2y = c2 * (a2x - a1x);

                        s1x = a0x + m1x;

                        a[aoffset + 0] = a0x;
                        a[aoffset + 1] = s1x - m2y;
                        a[aoffset + 2] = s1x + m2y;
                        return;
                    }

                    if(n == 4)
                    {
                        a0x = a[aoffset + 0];
                        a1x = a[aoffset + 1];
                        a2x = a[aoffset + 2];
                        a3x = a[aoffset + 3];

                        t1x = a0x + a2x;
                        t2x = a1x + a3x;

                        m2x = a0x - a2x;
                        m3y = a3x - a1x;

                        a[aoffset + 0] = t1x + t2x;
                        a[aoffset + 1] = m2x - m3y;
                        a[aoffset + 2] = t1x - t2x;
                        a[aoffset + 3] = m2x + m3y;
                        return;
                    }

                    if(n == 5)
                    {
                        offs = plan.plan[entryoffset + 7];

                        c1 = plan.precomputed[offs + 0];
                        c2 = plan.precomputed[offs + 1];
                        c3 = plan.precomputed[offs + 2];
                        c4 = plan.precomputed[offs + 3];
                        c5 = plan.precomputed[offs + 4];

                        t1x = a[aoffset + 1] + a[aoffset + 4];
                        t2x = a[aoffset + 2] + a[aoffset + 3];
                        t3x = a[aoffset + 1] - a[aoffset + 4];
                        t4x = a[aoffset + 3] - a[aoffset + 2];

                        t5x = t1x + t2x;
                        v0 = a[aoffset + 0] + t5x;
                        a[aoffset + 0] = v0;

                        m2x = c2 * (t1x - t2x);
                        m3y = c3 * (t3x + t4x);

                        s3y = m3y - c4 * t4x;
                        s5y = m3y + c5 * t3x;
                        s1x = v0 + c1 * t5x;

                        s2x = s1x + m2x;
                        s4x = s1x - m2x;

                        a[aoffset + 1] = s2x - s3y;
                        a[aoffset + 2] = s4x - s5y;
                        a[aoffset + 3] = s4x + s5y;
                        a[aoffset + 4] = s2x + s3y;
                        return;
                    }
                }

                if(plan.plan[entryoffset + 3] != fftbluesteinplan) return;
                //
                // Bluestein plan:
                // 1. multiply by precomputed coefficients
                // 2. make convolution: forward FFT, multiplication by precomputed FFT
                //    and backward FFT. backward FFT is represented as
                //
                //        invfft(x) = fft(x')'/M
                //
                //    for performance reasons reduction of inverse FFT to
                //    forward FFT is merged with multiplication of FFT components
                //    and last stage of Bluestein's transformation.
                // 3. post-multiplication by Bluestein factors
                //
                n = plan.plan[entryoffset + 1];
                var m = plan.plan[entryoffset + 4];
                offs = plan.plan[entryoffset + 7];

                for(i = stackptr + 2 * n; i <= stackptr + 2 * m - 1; i++)
                    plan.stackbuf[i] = 0;

                var offsp = offs + 2 * m;
                offsa = aoffset;
                var offsb = stackptr;

                double x;
                double y;
                double bx;
                double by;
                for(i = 0; i <= n - 1; i++)
                {
                    bx = plan.precomputed[offsp + 0];
                    by = plan.precomputed[offsp + 1];

                    x = a[offsa + 0];
                    y = a[offsa + 1];

                    plan.stackbuf[offsb + 0] = x * bx - y * -by;
                    plan.stackbuf[offsb + 1] = x * -by + y * bx;

                    offsp = offsp + 2;
                    offsa = offsa + 2;
                    offsb = offsb + 2;
                }

                ftbaseexecuteplanrec(ref plan.stackbuf, stackptr, plan, plan.plan[entryoffset + 5], stackptr + 2 * 2 * m);

                offsb = stackptr;
                offsp = offs;
                for(i = 0; i <= m - 1; i++)
                {
                    x = plan.stackbuf[offsb + 0];
                    y = plan.stackbuf[offsb + 1];

                    bx = plan.precomputed[offsp + 0];
                    by = plan.precomputed[offsp + 1];

                    plan.stackbuf[offsb + 0] = x * bx - y * by;
                    plan.stackbuf[offsb + 1] = -(x * by + y * bx);

                    offsb = offsb + 2;
                    offsp = offsp + 2;
                }

                ftbaseexecuteplanrec(ref plan.stackbuf, stackptr, plan, plan.plan[entryoffset + 5], stackptr + 2 * 2 * m);

                offsb = stackptr;
                offsp = offs + 2 * m;
                offsa = aoffset;

                for(i = 0; i <= n - 1; i++)
                {
                    x = plan.stackbuf[offsb + 0] / m;
                    y = -(plan.stackbuf[offsb + 1] / m);

                    bx = plan.precomputed[offsp + 0];
                    by = plan.precomputed[offsp + 1];

                    a[offsa + 0] = x * bx - y * -by;
                    a[offsa + 1] = x * -by + y * bx;

                    offsp = offsp + 2;
                    offsa = offsa + 2;
                    offsb = offsb + 2;
                }
                return;
            }

            /*************************************************************************
            Returns good factorization N=N1*N2.

            Usually N1<=N2 (but not always - small N's may be exception).
            if N1<>1 then N2<>1.

            Factorization is chosen depending on task type and codelets we have.

              -- ALGLIB --
                 Copyright 01.05.2009 by Bochkanov Sergey
            *************************************************************************/
            public static void ftbasefactorize(int n, out int n1, out int n2)
            {
                int j;

                n1 = 0;
                n2 = 0;

                //
                // try to find good codelet
                //
                if(n1 * n2 != n)
                    for(j = ftbasecodeletrecommended; j >= 2; j--)
                        if(n % j == 0)
                        {
                            n1 = j;
                            n2 = n / j;
                            break;
                        }

                //
                // try to factorize N
                //
                if(n1 * n2 != n)
                    for(j = ftbasecodeletrecommended + 1; j <= n - 1; j++)
                        if(n % j == 0)
                        {
                            n1 = j;
                            n2 = n / j;
                            break;
                        }

                //
                // looks like N is prime :(
                //
                if(n1 * n2 != n)
                {
                    n1 = 1;
                    n2 = n;
                }

                //
                // normalize
                //
                if(n2 != 1 || n1 == 1) return;
                n2 = n1;
                n1 = 1;
            }


            /*************************************************************************
            Is number smooth?

              -- ALGLIB --
                 Copyright 01.05.2009 by Bochkanov Sergey
            *************************************************************************/
            /*
                        public static bool ftbaseissmooth(int n)
                        {

                            for(var i = 2; i <= ftbasemaxsmoothfactor; i++)
                                while(n % i == 0)
                                    n /= i;
                            return n == 1;
                        }
            */


            /*************************************************************************
            Returns smallest smooth (divisible only by 2, 3, 5) number that is greater
            than or equal to max(N,2)

              -- ALGLIB --
                 Copyright 01.05.2009 by Bochkanov Sergey
            *************************************************************************/
            public static int ftbasefindsmooth(int n)
            {
                var best = 2;
                while(best < n) best = 2 * best;
                ftbasefindsmoothrec(n, 1, 2, ref best);
                return best;
            }


            /*************************************************************************
            Returns  smallest  smooth  (divisible only by 2, 3, 5) even number that is
            greater than or equal to max(N,2)

              -- ALGLIB --
                 Copyright 01.05.2009 by Bochkanov Sergey
            *************************************************************************/
            /*
                        public static int ftbasefindsmootheven(int n)
                        {
                            var best = 2;
                            while(best < n) best = 2 * best;
                            ftbasefindsmoothrec(n, 2, 2, ref best);
                            return best;
                        }
            */


            /*************************************************************************
            Returns estimate of FLOP count for the FFT.

            It is only an estimate based on operations count for the PERFECT FFT
            and relative inefficiency of the algorithm actually used.

            N should be power of 2, estimates are badly wrong for non-power-of-2 N's.

              -- ALGLIB --
                 Copyright 01.05.2009 by Bochkanov Sergey
            *************************************************************************/
            /*
                        public static double ftbasegetflopestimate(int n)
                        {
                            return ftbaseinefficiencyfactor * (4 * n * Math.Log(n) / Math.Log(2) - 6 * n + 8);
                        }
            */


            /*************************************************************************
            Recurrent subroutine for the FFTGeneratePlan:

            PARAMETERS:
                N                   plan size
                IsReal              whether input is real or not.
                                    subroutine MUST NOT ignore this flag because real
                                    inputs comes with non-initialized imaginary parts,
                                    so ignoring this flag will result in corrupted output
                HalfOut             whether full output or only half of it from 0 to
                                    floor(N/2) is needed. This flag may be ignored if
                                    doing so will simplify calculations
                Plan                plan array
                PlanSize            size of used part (in integers)
                PrecomputedSize     size of precomputed array allocated yet
                PlanArraySize       plan array size (actual)
                TmpMemSize          temporary memory required size
                BluesteinMemSize    temporary memory required size

              -- ALGLIB --
                 Copyright 01.05.2009 by Bochkanov Sergey
            *************************************************************************/
            private static void ftbasegenerateplanrec(int n, int tasktype, ftplan plan, ref int plansize, ref int precomputedsize, ref int planarraysize, ref int tmpmemsize, ref int stackmemsize, int stackptr)
            {
                //
                // prepare
                //
                if(plansize + ftbaseplanentrysize > planarraysize)
                    fftarrayresize(ref plan.plan, ref planarraysize, 8 * planarraysize);
                var entryoffset = plansize;
                const int esize = ftbaseplanentrysize;
                plansize = plansize + esize;

                //
                // if N=1, generate empty plan and exit
                //
                if(n == 1)
                {
                    plan.plan[entryoffset + 0] = esize;
                    plan.plan[entryoffset + 1] = -1;
                    plan.plan[entryoffset + 2] = -1;
                    plan.plan[entryoffset + 3] = fftemptyplan;
                    plan.plan[entryoffset + 4] = -1;
                    plan.plan[entryoffset + 5] = -1;
                    plan.plan[entryoffset + 6] = -1;
                    plan.plan[entryoffset + 7] = -1;
                    return;
                }

                //
                // generate plans
                //
                int n1;
                int n2;
                ftbasefactorize(n, out n1, out n2);
                if(tasktype == ftbasecffttask || tasktype == ftbaserffttask)
                {
                    //
                    // complex FFT plans
                    //
                    if(n1 != 1)
                    {
                        //
                        // Cooley-Tukey plan (real or complex)
                        //
                        // Note that child plans are COMPLEX
                        // (whether plan itself is complex or not).
                        //
                        tmpmemsize = Math.Max(tmpmemsize, 2 * n1 * n2);
                        plan.plan[entryoffset + 0] = esize;
                        plan.plan[entryoffset + 1] = n1;
                        plan.plan[entryoffset + 2] = n2;

                        plan.plan[entryoffset + 3] = tasktype == ftbasecffttask
                                                         ? fftcooleytukeyplan
                                                         : fftrealcooleytukeyplan;

                        plan.plan[entryoffset + 4] = 0;
                        plan.plan[entryoffset + 5] = plansize;
                        ftbasegenerateplanrec(n1, ftbasecffttask, plan, ref plansize, ref precomputedsize, ref planarraysize,
                                              ref tmpmemsize, ref stackmemsize, stackptr);
                        plan.plan[entryoffset + 6] = plansize;
                        ftbasegenerateplanrec(n2, ftbasecffttask, plan, ref plansize, ref precomputedsize, ref planarraysize,
                                              ref tmpmemsize, ref stackmemsize, stackptr);
                        plan.plan[entryoffset + 7] = -1;
                        return;
                    }

                    if(n == 2 || n == 3 || n == 4 || n == 5)
                    {
                        //
                        // hard-coded plan
                        //
                        plan.plan[entryoffset + 0] = esize;
                        plan.plan[entryoffset + 1] = n1;
                        plan.plan[entryoffset + 2] = n2;
                        plan.plan[entryoffset + 3] = fftcodeletplan;
                        plan.plan[entryoffset + 4] = 0;
                        plan.plan[entryoffset + 5] = -1;
                        plan.plan[entryoffset + 6] = -1;
                        plan.plan[entryoffset + 7] = precomputedsize;

                        if(n == 3) precomputedsize = precomputedsize + 2;

                        if(n == 5) precomputedsize = precomputedsize + 5;
                        return;
                    }
                    //
                    // Bluestein's plan
                    //
                    // Select such M that M>=2*N-1, M is composite, and M's
                    // factors are 2, 3, 5
                    //
                    var k = 2 * n2 - 1;
                    var m = ftbasefindsmooth(k);
                    tmpmemsize = Math.Max(tmpmemsize, 2 * m);
                    plan.plan[entryoffset + 0] = esize;
                    plan.plan[entryoffset + 1] = n2;
                    plan.plan[entryoffset + 2] = -1;
                    plan.plan[entryoffset + 3] = fftbluesteinplan;
                    plan.plan[entryoffset + 4] = m;
                    plan.plan[entryoffset + 5] = plansize;

                    stackptr = stackptr + 2 * 2 * m;
                    stackmemsize = Math.Max(stackmemsize, stackptr);

                    ftbasegenerateplanrec(m, ftbasecffttask, plan, ref plansize, ref precomputedsize,
                                          ref planarraysize, ref tmpmemsize, ref stackmemsize, stackptr);

                    //stackptr = stackptr - 2 * 2 * m;
                    plan.plan[entryoffset + 6] = -1;
                    plan.plan[entryoffset + 7] = precomputedsize;
                    precomputedsize = precomputedsize + 2 * m + 2 * n;
                    return;
                }
                if(tasktype != ftbaserfhttask) return;
                //
                // real FHT plans
                //
                if(n1 != 1)
                {
                    //
                    // Cooley-Tukey plan
                    //
                    //
                    tmpmemsize = Math.Max(tmpmemsize, 2 * n1 * n2);
                    plan.plan[entryoffset + 0] = esize;
                    plan.plan[entryoffset + 1] = n1;
                    plan.plan[entryoffset + 2] = n2;
                    plan.plan[entryoffset + 3] = fhtcooleytukeyplan;
                    plan.plan[entryoffset + 4] = 0;
                    plan.plan[entryoffset + 5] = plansize;

                    ftbasegenerateplanrec(n1, tasktype, plan, ref plansize, ref precomputedsize, ref planarraysize,
                                          ref tmpmemsize, ref stackmemsize, stackptr);

                    plan.plan[entryoffset + 6] = plansize;

                    ftbasegenerateplanrec(n2, tasktype, plan, ref plansize, ref precomputedsize, ref planarraysize,
                                          ref tmpmemsize, ref stackmemsize, stackptr);

                    plan.plan[entryoffset + 7] = -1;
                    return;
                }
                //
                // N2 plan
                //
                plan.plan[entryoffset + 0] = esize;
                plan.plan[entryoffset + 1] = n1;
                plan.plan[entryoffset + 2] = n2;
                plan.plan[entryoffset + 3] = fhtn2plan;
                plan.plan[entryoffset + 4] = 0;
                plan.plan[entryoffset + 5] = -1;
                plan.plan[entryoffset + 6] = -1;
                plan.plan[entryoffset + 7] = -1;

                if(n == 2 || n == 3 || n == 4 || n == 5)
                {
                    //
                    // hard-coded plan
                    //
                    plan.plan[entryoffset + 0] = esize;
                    plan.plan[entryoffset + 1] = n1;
                    plan.plan[entryoffset + 2] = n2;
                    plan.plan[entryoffset + 3] = fhtcodeletplan;
                    plan.plan[entryoffset + 4] = 0;
                    plan.plan[entryoffset + 5] = -1;
                    plan.plan[entryoffset + 6] = -1;
                    plan.plan[entryoffset + 7] = precomputedsize;

                    if(n == 3) precomputedsize = precomputedsize + 2;

                    if(n == 5) precomputedsize = precomputedsize + 5;
                    return;
                }
                return;
            }


            /*************************************************************************
            Recurrent subroutine for precomputing FFT plans

              -- ALGLIB --
                 Copyright 01.05.2009 by Bochkanov Sergey
            *************************************************************************/
            private static void ftbaseprecomputeplanrec(ftplan plan, int entryoffset, int stackptr)
            {

                if(plan.plan[entryoffset + 3] == fftcooleytukeyplan ||
                    plan.plan[entryoffset + 3] == fftrealcooleytukeyplan ||
                    plan.plan[entryoffset + 3] == fhtcooleytukeyplan)
                {
                    ftbaseprecomputeplanrec(plan, plan.plan[entryoffset + 5], stackptr);
                    ftbaseprecomputeplanrec(plan, plan.plan[entryoffset + 6], stackptr);
                    return;
                }

                int n;
                int offs;
                if(plan.plan[entryoffset + 3] == fftcodeletplan || plan.plan[entryoffset + 3] == fhtcodeletplan)
                {
                    var n1 = plan.plan[entryoffset + 1];
                    var n2 = plan.plan[entryoffset + 2];
                    n = n1 * n2;

                    const double pi23 = Consts.pi2 / 3;
                    if(n == 3)
                    {
                        offs = plan.plan[entryoffset + 7];
                        plan.precomputed[offs + 0] = Math.Cos(pi23) - 1;
                        plan.precomputed[offs + 1] = Math.Sin(pi23);
                        return;
                    }

                    if(n == 5)
                    {
                        offs = plan.plan[entryoffset + 7];
                        const double pi25 = Consts.pi2 / 5;
                        const double pi45 = pi25 * 2;
                        plan.precomputed[offs + 0] = (Math.Cos(pi25) + Math.Cos(pi45)) / 2 - 1;
                        plan.precomputed[offs + 1] = (Math.Cos(pi25) - Math.Cos(pi45)) / 2;
                        plan.precomputed[offs + 2] = -Math.Sin(pi25);
                        plan.precomputed[offs + 3] = -(Math.Sin(pi25) + Math.Sin(pi45));
                        plan.precomputed[offs + 4] = Math.Sin(pi25) - Math.Sin(pi45);
                        return;
                    }
                }

                if(plan.plan[entryoffset + 3] != fftbluesteinplan) return;

                ftbaseprecomputeplanrec(plan, plan.plan[entryoffset + 5], stackptr);

                n = plan.plan[entryoffset + 1];
                var m = plan.plan[entryoffset + 4];
                offs = plan.plan[entryoffset + 7];

                int i;
                for(i = 0; i < 2 * m; i++)
                    plan.precomputed[offs + i] = 0;

                var pin = Math.PI / n;
                for(i = 0; i < n; i++)
                {
                    var bx = Math.Cos(pin * i * i);
                    var by = Math.Sin(pin * i * i);
                    plan.precomputed[offs + 2 * i + 0] = bx;
                    plan.precomputed[offs + 2 * i + 1] = by;
                    plan.precomputed[offs + 2 * m + 2 * i + 0] = bx;
                    plan.precomputed[offs + 2 * m + 2 * i + 1] = by;

                    if(i <= 0) continue;
                    plan.precomputed[offs + 2 * (m - i) + 0] = bx;
                    plan.precomputed[offs + 2 * (m - i) + 1] = by;
                }

                ftbaseexecuteplanrec(ref plan.precomputed, offs, plan, plan.plan[entryoffset + 5], stackptr);
                return;
            }


            /*************************************************************************
            Twiddle factors calculation

              -- ALGLIB --
                 Copyright 01.05.2009 by Bochkanov Sergey
            *************************************************************************/
            private static void ffttwcalc(ref double[] a, int aoffset, int n1, int n2)
            {
                var n = n1 * n2;
                var pin = -Consts.pi2 / n;
                var twbasexm1 = Math.Sin(.5 * pin);
                twbasexm1 *= -2 * twbasexm1;
                //twbasexm1 = -(2 * math.sqr(Math.Sin(0.5 * v)));
                var twbasey = Math.Sin(pin);
                double twrowxm1 = 0;
                double twrowy = 0;
                for(var i = 0; i <= n2 - 1; i++)
                {
                    double twxm1 = 0;
                    double twy = 0;
                    double tmpy;
                    double tmpx;
                    int j;
                    for(j = 0; j <= n1 - 1; j++)
                    {
                        var idx = i * n1 + j;
                        var offs = aoffset + 2 * idx;
                        var x = a[offs + 0];
                        var y = a[offs + 1];
                        tmpx = x * twxm1 - y * twy;
                        tmpy = x * twy + y * twxm1;
                        a[offs + 0] = x + tmpx;
                        a[offs + 1] = y + tmpy;

                        //
                        // update Tw: Tw(new) = Tw(old)*TwRow
                        //
                        if(j >= n1 - 1) continue;
                        if(j % ftbaseupdatetw == 0)
                        {
                            pin = -(Consts.pi2 * i * (j + 1) / n);
                            twxm1 = Math.Sin(.5 * pin);
                            twxm1 *= -2 * twxm1;
                            //twxm1 = -(2 * math.sqr(Math.Sin(0.5 * v)));
                            twy = Math.Sin(pin);
                        }
                        else
                        {
                            tmpx = twrowxm1 + twxm1 * twrowxm1 - twy * twrowy;
                            tmpy = twrowy + twxm1 * twrowy + twy * twrowxm1;
                            twxm1 += tmpx;
                            twy += tmpy;
                        }
                    }

                    //
                    // update TwRow: TwRow(new) = TwRow(old)*TwBase
                    //
                    if(i >= n2 - 1) continue;
                    if(j % ftbaseupdatetw == 0)
                    {
                        pin = -(Consts.pi2 * (i + 1) / n);
                        twrowxm1 = Math.Sin(.5 * pin);
                        twrowxm1 *= -2 * twrowxm1;
                        //twrowxm1 = -(2 * math.sqr(Math.Sin(0.5 * v)));
                        twrowy = Math.Sin(pin);
                    }
                    else
                    {
                        tmpx = twbasexm1 + twrowxm1 * twbasexm1 - twrowy * twbasey;
                        tmpy = twbasey + twrowxm1 * twbasey + twrowy * twbasexm1;
                        twrowxm1 = twrowxm1 + tmpx;
                        twrowy = twrowy + tmpy;
                    }
                }
            }


            /*************************************************************************
            Linear transpose: transpose complex matrix stored in 1-dimensional array

              -- ALGLIB --
                 Copyright 01.05.2009 by Bochkanov Sergey
            *************************************************************************/
            private static void internalcomplexlintranspose(ref double[] a, int m, int n, int astart, ref double[] buf)
            {
                ffticltrec(ref a, astart, n, ref buf, 0, m, m, n);
                var lv_I1 = 0 - astart;
                for(var i = astart; i <= astart + 2 * m * n - 1; i++)
                    a[i] = buf[i + lv_I1];
            }


            /*************************************************************************
            Linear transpose: transpose real matrix stored in 1-dimensional array

              -- ALGLIB --
                 Copyright 01.05.2009 by Bochkanov Sergey
            *************************************************************************/
            private static void internalreallintranspose(ref double[] a, int m, int n, int astart, ref double[] buf)
            {

                fftirltrec(ref a, astart, n, ref buf, 0, m, m, n);
                var i1 = 0 - astart;
                for(var i = astart; i <= astart + m * n - 1; i++)
                    a[i] = buf[i + i1];
            }


            /*************************************************************************
            Recurrent subroutine for a InternalComplexLinTranspose

            Write A^T to B, where:
            * A is m*n complex matrix stored in array A as pairs of real/image values,
              beginning from AStart position, with AStride stride
            * B is n*m complex matrix stored in array B as pairs of real/image values,
              beginning from BStart position, with BStride stride
            stride is measured in complex numbers, i.e. in real/image pairs.

              -- ALGLIB --
                 Copyright 01.05.2009 by Bochkanov Sergey
            *************************************************************************/
            private static void ffticltrec(ref double[] a, int astart, int astride, ref double[] b,
                                           int bstart, int bstride, int m, int n)
            {
                if(m == 0 || n == 0) return;
                if(Math.Max(m, n) <= 8)
                {
                    var m2 = 2 * bstride;
                    for(var i = 0; i <= m - 1; i++)
                    {
                        var idx1 = bstart + 2 * i;
                        var idx2 = astart + 2 * i * astride;
                        for(var j = 0; j <= n - 1; j++)
                        {
                            b[idx1 + 0] = a[idx2 + 0];
                            b[idx1 + 1] = a[idx2 + 1];
                            idx1 = idx1 + m2;
                            idx2 = idx2 + 2;
                        }
                    }
                    return;
                }
                if(n > m)
                {
                    //
                    // New partition:
                    //
                    // "A^T -> B" becomes "(A1 A2)^T -> ( B1 )
                    //                                  ( B2 )
                    //
                    var n1 = n / 2;
                    if(n - n1 >= 8 && n1 % 8 != 0)
                        n1 += 8 - n1 % 8;
                    Contract.Assert(n - n1 > 0);
                    ffticltrec(ref a, astart, astride, ref b, bstart, bstride, m, n1);
                    ffticltrec(ref a, astart + 2 * n1, astride, ref b, bstart + 2 * n1 * bstride, bstride, m, n - n1);
                }
                else
                {
                    //
                    // New partition:
                    //
                    // "A^T -> B" becomes "( A1 )^T -> ( B1 B2 )
                    //                     ( A2 )
                    //
                    var m1 = m / 2;
                    if(m - m1 >= 8 && m1 % 8 != 0)
                        m1 += 8 - m1 % 8;
                    Contract.Assert(m - m1 > 0);
                    ffticltrec(ref a, astart, astride, ref b, bstart, bstride, m1, n);
                    ffticltrec(ref a, astart + 2 * m1 * astride, astride, ref b, bstart + 2 * m1, bstride, m - m1, n);
                }
            }


            /*************************************************************************
            Recurrent subroutine for a InternalRealLinTranspose


              -- ALGLIB --
                 Copyright 01.05.2009 by Bochkanov Sergey
            *************************************************************************/
            private static void fftirltrec(ref double[] a, int astart, int astride, ref double[] b, int bstart,
                                           int bstride, int m, int n)
            {
                if(m == 0 || n == 0) return;

                if(Math.Max(m, n) <= 8)
                {
                    for(var i = 0; i <= m - 1; i++)
                    {
                        var idx1 = bstart + i;
                        var idx2 = astart + i * astride;
                        for(var j = 0; j <= n - 1; j++)
                        {
                            b[idx1] = a[idx2];
                            idx1 = idx1 + bstride;
                            idx2 = idx2 + 1;
                        }
                    }
                    return;
                }
                if(n > m)
                {
                    //
                    // New partition:
                    //
                    // "A^T -> B" becomes "(A1 A2)^T -> ( B1 )
                    //                                  ( B2 )
                    //
                    var n1 = n / 2;
                    if(n - n1 >= 8 && n1 % 8 != 0)
                        n1 += 8 - n1 % 8;
                    Contract.Assert(n - n1 > 0);
                    fftirltrec(ref a, astart, astride, ref b, bstart, bstride, m, n1);
                    fftirltrec(ref a, astart + n1, astride, ref b, bstart + n1 * bstride, bstride, m, n - n1);
                }
                else
                {
                    //
                    // New partition:
                    //
                    // "A^T -> B" becomes "( A1 )^T -> ( B1 B2 )
                    //                     ( A2 )
                    //
                    var m1 = m / 2;
                    if(m - m1 >= 8 && m1 % 8 != 0)
                        m1 += 8 - m1 % 8;
                    Contract.Assert(m - m1 > 0);
                    fftirltrec(ref a, astart, astride, ref b, bstart, bstride, m1, n);
                    fftirltrec(ref a, astart + m1 * astride, astride, ref b, bstart + m1, bstride, m - m1, n);
                }
            }


            /*************************************************************************
            recurrent subroutine for FFTFindSmoothRec

              -- ALGLIB --
                 Copyright 01.05.2009 by Bochkanov Sergey
            *************************************************************************/
            private static void ftbasefindsmoothrec(int n, int seed, int leastfactor, ref int best)
            {
                Contract.Assert(ftbasemaxsmoothfactor <= 5, "FTBaseFindSmoothRec: internal error!");
                if(seed >= n)
                {
                    best = Math.Min(best, seed);
                    return;
                }
                if(leastfactor <= 2)
                    ftbasefindsmoothrec(n, seed * 2, 2, ref best);
                if(leastfactor <= 3)
                    ftbasefindsmoothrec(n, seed * 3, 3, ref best);
                if(leastfactor <= 5)
                    ftbasefindsmoothrec(n, seed * 5, 5, ref best);
            }


            /*************************************************************************
            Internal subroutine: array resize

              -- ALGLIB --
                 Copyright 01.05.2009 by Bochkanov Sergey
            *************************************************************************/
            private static void fftarrayresize(ref int[] a, ref int asize, int newasize)
            {
                var tmp = new int[asize];
                for(var i = 0; i <= asize - 1; i++)
                    tmp[i] = a[i];

                a = new int[newasize];
                for(var i = 0; i <= asize - 1; i++)
                    a[i] = tmp[i];

                asize = newasize;
            }


            /*************************************************************************
            Reference FHT stub
            *************************************************************************/
            private static void reffht(ref double[] a, int n, int offs)
            {
                Contract.Assert(n > 0, "RefFHTR1D: incorrect N!");
                var buf = new double[n];
                var pin = Consts.pi2 / n;
                for(var i = 0; i < n; i++)
                {
                    double v = 0;
                    int j;
                    for(j = 0; j <= n - 1; j++)
                    {
                        var ij = i * j;
                        v += a[offs + j] * (Math.Cos(pin * ij) + Math.Sin(pin * ij));
                    }
                    buf[i] = v;
                }

                for(var i = 0; i <= n - 1; i++)
                    a[offs + i] = buf[i];
            }

            #region Nested type: ftplan

            public class ftplan
            {
                public int[] plan;
                public double[] precomputed;
                public double[] stackbuf;
                public double[] tmpbuf;

                public ftplan()
                {
                    plan = new int[0];
                    precomputed = new double[0];
                    tmpbuf = new double[0];
                    stackbuf = new double[0];
                }
            } ;

            #endregion
        }

    }

}